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ON  GENERATION  OF  RANDOM  NUMBERS  WITH  SPECIFIED 
DISTRIBUTIONS  OR  DENSITIES 

INTRODUCTION 

The  generation  of  random  numbers  with  a  specified  probability  density 
function  or  a  specified  cumulative  distribution  function  is  a  frequent 
occurrence  in  the  simulation  of  signal  processing  techniques  that  are 
difficult  or  impossible  to  evaluate  analytically.  Accordingly,  it  is  of 
interest  to  be  able  to  generate  easily  and  efficiently  such  random  numbers. 
It  is  assumed  that  a  uniform  random  number  generator  is  already  available; 
that  is,  independent  random  variables  with  a  probability  density  function 
equal  to  unity  over  the  range  (0,1)  can  be  generated.  Validation  of  this 
assumption  for  a  particular  random  number  generator  must  be  confirmed  by 
numerical  tests,  several  of  which  are  tested  herein.  This  report  is  an 
extension  and  elaboration  of  [1]. 
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NONLINEAR  DISTORTION  OF  A  SINGLE  UNIFORM  RANDOM  VARIABLE 

Suppose  that  random  variable  x  is  uniformly  distributed  on  (0,1).  Let  the 
nonlinear  distortion  of  x  yield  the  random  variable  y  given  by 


pMx> 

y  =  g(x)  =J  or 

K<x> 

where 


9 -i (  )  is  a  monotonical ly  increasing  function  of  its  argument, 

9d(  )  is  a  monotonical ly  decreasing  function  of  its  argument.  (2^ 

This  restriction  of  the  nonlinear  distortion  to  be  monotonic  (either 
increasing  or  decreasing)  is  not  necessary;  it  is  done  here  only  for 
simplicity.  Then  since  random  variable  x  is  uniformly  distributed  on  (0,1), 
for  a  fixed  (nonrandom)  t  in  the  range  (0,1), 


t  =  Prob 


{*<  tj . 


y  <  9i(t) 


=  Prob  <  or 


y  >  gd(t) 


g^x)  <  gi  (t) 


9d(>0  >  gdU) 


Py(9-j  (t)  ) 


1  -  Py(9d(t)) 


where  we  used  (2),  (1),  and  denoted  the  cumulative  distribution  function  of 
random  variable  y  by  Py(  ).  we  rewrite  (3)  as 

t  =  P  y  ( 9  i  ( t ) ) 

or  (4 

l-t=P(gd(t))  , 


and  observe  that  cumulative  distribution  function  Py(  )  is  a  monotonical ly 
increasing  function  of  its  argument. 
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Now,  for  arbitrary  monotonic  function  h(  ),  we  denote  its  inverse  function 
by  h(  );  that  is,* 

h(h(u) )  =  u,  h(fi(v) )  =  v  (5) 

Applying  monotonic  function  P  (  )  to  both  sides  of  (4),  and  using  (5),  there 
follows 


g^t)  =  py(t) 
or 

gd(t)  *  fy1  "  * 

Employing  (6)  in  (1)  final ly  yields  the  random  variable 


y  =  9(x)  = 


fow -fy(X>  " 

1  or 

[gd(x)  =  Py(i  -  x)  j 


(6) 


(7) 


as  the  desired  result  of  this  nonlinear  transformation  of  random  variable  x. 
If  cumulative  distribution  function  P  (  )  of  random  variable  y  is  specified 
and  desired,  and  if  random  variable  x  is  uniformly  distributee  on  (0,1), 

(7)  tells  us  that  we  must  evaluate  the  inverse  of  the  given  cumulative 
distribution  function  and  use  it  as  the  nonlinear  transformation  of  either 
random  variable  x  or  random  variable  1-x,  depending  on  whether  we  want  a 
monotonical ly  increasing  or  monotonically  decreasing  transformation, 
respectively.  The  key  element  to  this  approach  is  the  ease  with  which  the 
inverse  function  T*  (  )  can  be  computed. 


*  Some  additional  relationships  among  functions  and  their  inverses  are 
presented  in  appendix  A. 
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EXAMPLES  OF  SINGLE-VARIATE  DISTORTION 


The  following  examples  have  been  adjusted  to  a  convenient  scale,  such  as 
zero  mean  or  unit  variance.  By  addition  of  a  constant  and/or  multiplication 
by  a  scale  factor,  alternative  desired  ranges  can  be  realized.  The  times  of 
execution  given  below  were  obtained  for  the  Hewlett-Packard  9845B  Desk 
Calculator  equipped  with  the  Fast  Processor  Upgrade  Kit;  they  include  the  time 
required  to  generate  the  uniform  random  variable  x.  Results  were  obtained  by 
averaging  over  1000  independent  trials.  Loop  counters  were  declared  INTEGER 
for  maximum  speed. 

EXPONENTIAL  DENSITY 

The  desired  probability  density  function  of  random  variable  y  is 
exponential*: 

py(u)  =  exp(-u)  for  u  >  0  (8) 

The  corresponding  cumulative  distribution  function  is 

t 

P  (t)  =  f  du  p  (u)  =  1  -  exp(-t)  for  t  >  0  .  (9) 

The  characteristic  function  of  random  variable  y  is  (for  future  reference) 

+•  i 

fy  (?)  =  J  du  exp(iju)  py(u)  =  (1  -  i?)  .  (10) 

-  OO 

In  order  to  find  the  inverse  function  to  (9),  we  let 

u  =  Py(t)  =  1  -  exp(-t)  .  (11) 

Solving  the  left  relation  for  t,  and  then  solving  the  outside  relation  for  t, 
we  get,  respectively, 

*  The  probability  density  functiorwand  cumulative  distribution  functions  are 
zero  in  the  unspecified  regions;  for  example,  Py(u)  =  0  for  u  <  0. 
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2 

Py(u)  =  u  exp(-u  12)  for  u  >  0  , 

Py(u)  =  1  -  exp(-u^/2)  for  u  >  0  .  (15) 

Via  a  procedure  similar  to  (11)  and  (12),  there  follows 


Py(u)  «  y-2  in(l-u) 


0  <  u  < 


The  nonlinear  transformation  that  yields  Rayleigh  random  variables  is 


The  time  of  execution  of  (17)  is  2.5  msec. 


(16) 


(17) 


5 


Gaussian  density 


Py( u )  =  (2ir)  1/2  exp(-u2/2)  for  all  u,  \ 

py(u)  =  jf  dt  (2*r1/2  exp(-t2/2)  s  $(u)  ,  (18)' 

-o 0 


fy(f)  =  exp  (-fl 2)  . 

The  inverse  function  is  [2;  26.2,  26.2.22,  26.2.23] 

Py(u)  *  f(u)  *  -f(l-u)  for  0  <  u  <  1  .  (19) 


The  nonlinear  transformation  that  yields  Gaussian  random  variables  from 
uniformly  distributed  ones  is  then,  from  (7)  and  (19), 


The  time  of  execution  of  (20)  is  13.2  msec. 


(20) 


A  much  better  approach,  for  this  particular  case  of  generation  of  Gaussian 
random  variables,  is  as  follows:  let  x^  and  ^  &e  two  independent  random 
variables,  uniformly  distributed  on  (0,1).  Then  according  to  the  previous 
example,  the  two  independent  random  variables 

r  in  x^  ,  9  =  2*  x2  ,  (21) 

have,  respectively,  the  probability  density  functions 


2 

Pr(u)  =  u  exp(-u  /2)  for  u  >  0  , 

p  (u)  =  (2*)"1  for  0  <  u  <  2w  .  (22) 
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Now  define  two  new  random  variables  by  the  nonlinear  transformations 


y1  =  r  cos  o,  =  r  sin 


(23) 


The  joint  characteristic  function  of  y1  and  y2  is  [2;  y.1.21  and  11.4.29] 

E^expOfyj  +  it»y2j]  =  E[exp(i$r  cos  o  +  iwr  sin  ejj 
►  s 

du  dv  Pr(u)  ,pQ(v)  exp  (iju  cos  v  +  iyu  sin  v) 


-JT' 

-o© 

4-flO 

*  l 


du  u  exp(-u*72) 


dv 


27  exp(iu  (f  cos  v  +  u  sin  v)) 


■w 

£  du  u  exp(-u2/2)  J 


- 1  e2  1  2 v  I  \ 

exP(-  J  5  ~  2  v  ' 


OF77) 


(24) 


siaiXr 


andom  variables,  each  with 


Thus  yj  and  y2  are  independent  Gaus 

zero-mean  and  unit  variance.  The  time  of.  execution  of  (21)  and  (23)  is 
5.4  msec  per  random  variable  (actually  loV  msec  for  a  pair  of  independent 
Gaussian  random  variables.)  This  5.4  msec^s  considerably  less  than  the 
i3.2  msec  required  of  (20),  which  also  requires  a  special  function 
definition.  A  more  general  distortion  than  1 )— (2 3 )  is  considered  in 

appendix  B.  \ 

\ 

Another  alternative  for  the  generation  of  approximately  Gaussian  random 
variables  is  to  sum  N  independent  random  variable^,  uniformly  distributed  over 
(0,1).  By  subtraction  of  the  constant  N/2  and  scal'ing  by  ^12/n;  a  zero-mean 
unit-variance  random  variable  can  be  generated  which,  however,  is  limited  to 
the  finite  range  (-V3N1,  ).  A  table  of  execution  times  and  range  values  is 

given  below. 
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Table  1.  Execution  Times  for  Sum  of  N  Independent  Random  Variables 


Execution  Time  (msec) 


Range  of  Values 


4  \ 

2.8 

-3.46, 

3.46 

5 

\ 

3.3 

-3 .87, 

3.87 

6 

3.9 

-4.24, 

4.24 

7 

\ 

\ 

4.4 

-4.58, 

4.58 

8 

\  4.9 

-4.90, 

4.90 

9 

5v4 

-5.20, 

5.20 

10 

5.9\ 

-5.48, 

5.48 

20 

11.0 

-7.75, 

7.75 

30 

16.2 

-9.49, 

9.49 

The  characteristic  function  of  this  random  variabY&is  [sin(Rf)/(Rf)]N  where 
R  =  V3/N' ,  and  the  normalized  fourth-cumulant  of  the  random  variable  generated 
by  this  summation  procedure  is  -1.2/N.  If  the  non-Gaussi^anness  can  be 
tolerated,  this  summation  procedure  is  then  a  viable  alternative  to  (21)  and 
(23)  in  terms  of  execution  time,  for  N  <  9. 


CAUCHY  DENSITY 


P»  =  7 - -—7  for  all  u, 

y  *  i  +  / 


py(u)  -  \  +  7  arc  tan(u) , 


fy(?)  =  exp  (-|f|)  .  (25) 

This  is  the  probability  density  function  of  the  ratio  of  two  zero-mean 
independent  Gaussian  random  variables.  The  inverse  function  to  the  cumulative 
distribution  function  is 


£y(u)  =  tan 


[*  f  •  ?)] 


for  0  <  u  <  1 


(26) 


and  the  nonlinear  transformation  is 
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y  = 


(27) 


The  execution  time  of  (27)  is  2.8  msec.  Equivalently,  the  transformation 
tan(irx)  would  realize  the  desired  probability  density  function,  although  it  is 
not  monotonic  in  x  over  (0,1). 


RECTIFIED  CAUCHY  DENSITY 


for  u  >  0  , 

for  u  >  0  , 

for  0  <  u  <  1  , 


(28) 


\ 


The  execution  time  of  (28)  is  2.7  msec. 
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REALIZATION  OF  DISTRIBUTION  VIA  COMBINATION  OF  SEVERAL  RANDOM  VARIABLES 

If  several  independent  random  variables  are  added,  their  characteristic 
functions  multiply.  So  if  a  specified  probability  density  function  or 
cumulative  distribution  function  has  a  characteristic  function  which  can  be 
broken  down  into  easily  realizable  terms,  this  observation  can  be  utilized  to 
generate  complicated  distributions  with  relative  ease.  For  example,  consider 
the  following. 

CHI-SQUARED  DENSITY  WITH  2N  DEGREES  OF  FREEDOM 

This  variate  is  normally  thought  of  as  being  generated  by  summing  the 
squares  of  2N  zero-mean  unit-variance  independent  Gaussian  random  variables. 
One  half  of  this  sum  has  the  probability  density  function 


The  inverse  function,  Py(  )t  to  (30)  is  not  available  in  closed  form  for 
N  >  2,  so  that  recourse  to  (7)  is  not  reasonable. 


However,  the  characteristic  function  corresponding  to  (29)  is 


fy(f)  =  (1  -  i?fN  •  (31) 

But  this  is  (10)  multiplied  by  itself  N  times.  Then  (14)  reveals  that  we  can 
generate  y  according  to  a  sum  of  N  independent  variates: 


The  latter  form  in  (32)  is  preferable  computationally;  it  involves  N-l 
multiplies  and  only  one  logarithm,  and  is  obviously  much  quicker  than 
10 
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generating  2N  Gaussian  random  variables  and  summing  their  squares.  The  time 
of  execution  of  (32)  is  given  below. 


Table  2.  Execution  Times  for  Chi-Squared  Variate  of  2N  Degrees  of  Freedom 


Execution  Time  (msec) 


If  we  were  to  use  the  nonlinear  distortion  procedure  in  (20)  for  generating 
Gaussian  random  variables,  and  square  and  add  2N  such  numbers  to  generate  a 
Chi-squared  variate,  the  time  of  execution  would  be  2N  (13.2)  msec.  For  the 
examples  in  table  2,  these  times  are  264,  528,  and  792  msec  respectively, 
which  are  far  greater  than  the  execution  times  for  (32). 


A  closely  related  random  variable  to  (32)  is 


F  r  N 

=  *  ~2^n]TTxn  *9^ 

L  ln-1  J. 


The  probability  density  function  of  random  variable  r  is 


Pr(u)  =  py(9(u))  jjjflu)  =  Py(u2/2)  u 


4.2N-exp(V/2)  for  u>Q 
2N~1(N-1).' 


The  cumulative  distribution  function  of  r  is 


Pr(u)  ■  1  -  exp(-u*72)  ■  for  u  >  0 

n=0 
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The  time  of  execution  of  (33)  is  given  in  table  3. 

Table  3.  Execution  Times  for  Square  Root  of  Chi -Squared  Variate 


Another  example  of  this  combination  of  variates  is  furnished  by  the 
fol lowing. 

Q  DISTRIBUTION  AND  RICE  DENSITY 

Suppose  that  yj  and  y2  are  tw0  independent  zero-mean  unit-variance 
Gaussian  random  variables,  as  generated  via  (21)  and  (23).  Then  for  a 
constant  d,  the  random  variable 


2  -  I  [<yl  *  d)2  *  y2]“  ?  I?  *  2dyl  *  A  *  y0 

=  y  +  2dr  cos  o  +  r^J 


2dr  cos  o  +  r 
has  characteristic  function  [2;  9.1.21  and  11.4.29] 

,2 


(36) 


f,(5)  *  E{exp(i 


ifzj  =  Ejexp^if^ 


jj—  +■  dr  cos  e  + 

2ir 


m 


=  J  du  u  exp(-u^/2)  expjj^-  (d^  +  u^)J  ^  dv  ^  exp(i  J  d  u  cos  v) 

+<*  r  2  ”1 

«  exp( if d2/ 2)  J  du  u  exp  L  j-(l-if)J  JQ(f  d  u) 

Q  L  J 

-  d-'f'1  «p(#^fr)  • 


12 


N 

Execution  Time  (msec) 

10 

8.6 

r 

20 

14.7 

30 

20.8 

r 

r- 


» < 
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The  corresponding  probability  density  function  is 

ft 

• 

pz(u)  =  exp^  y~  -  u)  I0(dV25»)  for  u  >  0  , 

(38) 

>  fc  . 

as  may  be  confirmed  by  direct  Fourier  transformation  (according 

to  (10))  of 

ujill 

(38);  see  [2;  11.4.29].  The  cumulative  distribution  function  of  z  is 

Pz(u)  =  1  -  Q(d,V5TT  )  for  u  >  0  , 

(39) 

y 

where  Marcum's  Q  function  is  defined  as 

+«•  /  2  2\ 

Q(a,b)  -  f  dx  x  exp  2  )  IQ(ax) 

(40) 

r 

The  time  of  execution  of  (36)  and  (21)  is  7.2  msec;  there  is  no 
(23). 

need  to  employ 

ft 

A  closely  related  random  variable  to  (36)  is 

y  =  Y 2T  »  (d2  +  2dr  cos  e  +  r2^  . 

(41) 

ft 

The  probability  density  function  of  y  is  the  Rice  density 

Py(u)  =  Pz(u2/2)  u 

'  • 

r ; 

=  u  exp^  - — IQ(du)  for  u  >  0  » 

(42) 

and  the  cumulative  distribution  function  of  y  is 

• 

►  - 

P  (u)  =  1  -  Q(d,u)  for  u  >  0 

(43) 

The  time  of  execution  of  (41)  and  (21)  is  7.5  msec;  there  is  no 

need  to  employ 
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These  two  examples,  (36)  and  (41),  would  again  be  analytically  very 
difficult  to  realize  by  use  of  (7),  because  the  inverse  functions  to 
cumulative  distribution  functions  (39)  and  (43)  are  not  available  in  closed 
form.  Furthermore,  different  values  of  d  are  easily  accommodated  in  (36)  and 
(41),  whereas  the  inverses  to  (39)  and  (43)  would  necessarily  involve  d  as  a 
parameter. 

A  more  general  situation  is  encountered  for  the  following. 

Qm  DISTRIBUTION 

This  example  is  a  combination  of  the  last  two.  Let,  as  in  (36)  and  (21), 
v  -  j  [(yx  +  d)2  ♦  y^j  =  -|[d2  +  2dr  cos  e  +  r2]  ,  (44) 

where 

r  =^-2  in  xf  ,  e  =  2*  x2  .  (45) 

And  as  in  (32),  let 


M+l 

w  *  i  IT x 


m=3 


m 


,M+1 


(46) 


All  the  random  variables,  £xmj^  »  are  independent  and  uniformly  distributed 
on  (0,1).  Then  let  random  variable  z  be  the  sum  of  the  above  two: 


l  ?  2  f»M  1 

z  =  v+w  =  £[d  +  2dr  cos  ©  +  r  ]  -  in  <  TT  xm  j 

The  characteristic  function  of  random  variable  z  follows  upon  use  of 
(37)  and  (31): 


(47) 


j 


.  •'  1 
J 


-  i 
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fz(f)  =  d-if)41  exp(r  ir^)  •  (^) 

The  corresponding  probability  density  function  and  cumulative  distribution 
function  for  z  are 


The  execution  time  of  (45)  and  (47)  is  given  in  table  4  for  several  values  of 
M.  The  use  of  (7)  would  have  required  the  inverse  of  (49),  a  rather 
formidable  analytical  task  involving  the  two  parameters  d  and  M. 


Table  4.  Execution  Times  for  QM  Variate 


M 

Execution  Time  (msec) 

10 

14.6 

20 

20.6 

30 

26.7 

A  closely  related  random  variable  to  (47)  is 


(51) 


15 


There  follows 


py(u)  =  exp^  d  -g~U  j  IM_1(du)  for  u  >  0  , 

py(u)  =  1  -  QM(d,u)  for  u  >  0  .  (52) 

The  execution  times  for  (45)  and  (51)  were  0.4  msec  larger  than  those  given  in 
table  4. 


TR  6843 


PRODUCTS  OF  RANDOM  VARIABLES 

We  can  now  take  products  of  some  of  the  random  variables  above  and 
generate  additional  cases  of  complicated  probability  density  functions.  For 
example,  if  we  take  the  product  of  two  random  variables  as  given  by  (32), 

f  N  i  (N2 

Z  =  yxy2  =  Jin  J  J|~  xjj1*  >in  J|"  x^2)  ,  (53) 

|_n=l  J  n=l 

where  |x^jand  jxj^j  are  uniformly  distributed,  the  probability  density 
function  of  the  product  random  variable  is,  from  (29)  and  [4;  3.471  9j, 


Here  Ky(  )  iS  a  modified  Bessel  function  of  the  second  kind  and  order  v 

[2;  section  9.6];  it  is  even  in  v  [2;  9.6.6].  Execution  times  for  (53)  can  be 

determined  from  (32)  and  table  2. 

The  random  variable  w  =  2z^^,  where  z  is  given  by  (53),  has  the 
probability  density  function 

N1+N2-l 

pw(u)  =  (Nj-il:  (n2-i) :  kn1-n2(u)  for  u  >  0  *  (55) 

The  product  of  two  independent  zero-mean  Gaussian  random  variables,  as 
given  by  (21)  and  (23),  is  given  more  simply  by 
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2  12 

2  =  yly2  =  r  sin  9  cos  ©  =  j  r  sin  (2©)  =  — J^i(x^)  sin(4irx2)  .  (56) 

The  probability  density  function  of  z  can  be  determined  exactly  as  in  (54), 
wi th  the  result 

pz(u)  =  7  KQ(  [u| )  for  all  u  .  (57) 

The  4 it  sweep  of  the  sin  argument  in  (56)  is  unnecessary;  the  following  will 
accomplish  the  same  probability  density  function: 

z  =  jjn(xx)  cos(nx2)  ,  (58) 

where  x^  anc|  X£  are  uniform  over  (0,1).  The  time  of  execution  of  (58)  is 
6.0  msec. 
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GENERATION  FROM  A  NON  UNIFORMLY  DISTRIBUTED  RANDOM  VARIABLE 

Suppose  that  it  is  possible  to  generate  a  random  variable  x  with 
cumulative  distribution  function  Px(  ) ,  and  that  the  cumulative  distribution 
function  of 


y  =  9(x)  (5S 

is  desired  to  be  P  (  ).  Then  for  g  =  g-j,  a  monotonical ly  increasing 
function,  the  cumulative  distribution  function  of  x  is 

Px(u)  =  Prob  {x  <  u}  =  Prob  (g^x)  <  g^u)} 

-  Prob  fy  <  gi(u)J  •  Py(gi(u))  .  (6C 

Applying  inverse  function  P  to  both  sides  of  (60),  the  desired  nonlinearity 


9i(u)  =  Py(px(u))  •  0>U 

Equation  (7)  is  obviously  the  special  case  of  this  when  x  is  uniformly 
distributed.  Alternative  expressions  to  (61)  are  available  by  use  of  some  of 
the  relations  in  appendix  A. 

If  we  wart  to  use  distortion  g  =  gd>  a  monotonical ly  decreasing 
function,  to  generate  y,  we  have 


Px(u)  =  Prob  {x  <  u}  =  Prob  {gd(x)  >  gd(u)} 
-  Prob  [y  >  gd(u)}  =  1  -  PJgd(u)) 


Inverting  this  relation. 


9d(u>  -V1-px<u!> 


w 
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•  TESTS  OF  UNIFORM  RANDOM  NUMBER  GENERATOR 

A  critical  component  of  the  procedures  described  above  is  the  generator  of 
random  variables,  {xn}s  with  a  first-order  probability  density  function  that 
is  flat  over  (0,1), and  with  statistically  independent  samples.  Here  we  will 
investigate  several  tests  of  the  random  number  generator,  RND,  of  the 
Hewlett-Packard  9845  Desk  Calculator.  Inherent  in  this  investigation  is  the 
need  to  state  the  confidence  associated  with  a  particular  measurement  or 
estimate;  for  example,  see  [5]. 

CORRELATION  TEST 

The  sample  correlation  of  a  set  of  N  measurements  fyn3 1  is  defined 
here  as 


Rk  "  Nik  Vn-k  for  0  <  k  <  N  • 


We  presume  that  the  {y^  are  all  independent  with  a  flat  probability  density 
function 

Pv(u)  =  for  |U|  <VT  .  (6b 

y  2V? 

These  random  variables  can  be  obtained  from  uniform  {xn}  according  to  the 
linear  transformation 


-  *#(?«  -  0  : 


they  have  the  convenient  normalized  values 


E(yn)  =  0  ,  E(y‘)  =  1  . 


More  generally,  from  (66), 
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Under  these  assumptions,  it  is  readily  shown,  by  use  of  (69),  that  Rk  is 
unbiased,  that  is 


1  for  k  =  0 
0  otherwise 


E(Rk)  - 

and  that  the  standard  deviation  of  Rk  is 


(70) 


)(1.25N)"1/2 

o(R|/)  *  ]  -1/2 

k  '(N-k) 1/2 


for 

for 


k  =  0 
1  <  k  <  N 


1 


(71) 


Thus  we  expect  Rk  to  lie  within  2  or  3  standard  deviations,  a(Rk),  of 
E(Rk)  most  of  the  time;  that  is,  the  normalized  random  variable 

Rk  -  E(RJ 
rk  5  -50^ - 


(72) 


should  be  between  ±  2  most  of  the  time,  with  rare  excursions  to  ±  3,  if  the 
{xn}  are  truly  independent  [5].  In  table  5,  the  results  of  sample  runs  for 
N=100,  1000,  and  10000  are  listed.  They  furnish  no  reason  for  rejecting  the 

hypothesis  that  [xn|  are  independent.  Runs  for  other  sets  of  random  numbers 
yielded  results  very  similar  to  table  5.  An  alternative  test  on  the  whiteness 
of  {xn}  is  gi  ven  in  appendix  C. 

Table  5.  Sample  Correlation  Results  for  (72) 


Delay  k 

rk  for  N  «  100 

r^  for  N  =  1000 

r^  for  N  = 

0 

1.64 

0.78 

0.65 

1 

-1.04 

-0.91 

0.25 

2 

-0.26 

-0.68 

-0.02 

3 

-0.10 

1.39 

-0.13 

4 

1.44 

0.31 

0.27 

5 

-2.31 

-0.30 

0.47 

6 

-0.49 

-1.20 

-0.20 

7 

0.83 

0.81 

-1.40 

8 

0.17 

0.61 

-1.47 

9 

-0.61 

-0.77 

-0.91 

10 

-1.43 

-0.45 

0.48 
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MOMENT  TEST 

The  sample  moment  of  order  k  for  a  set  of  N  measurements  {yn}is 


The  mean  and  variance  of  uk  are 

E<l*k)  -  E(*nk) ,  var(„k)  -  i  var(ynk)  , 

and  can  be  obtained  from  (66)  and  (69).  The  random  variable 


(73) 


(74) 


-  E(u|>) 

mk  5  Std.  Dev.Tii'k')  (75) 

should  therefore  lie  mostly  in  the  range  ±2.  A  sample  result  is  given  in 
table  6,  which  again  confirms  the  *lat  probability  density  function  hypothesis 
in  (66).  (Althou^i  these  particular  runs  all  resulted  in  positive  numbers, 
other  sample  runs  resulted  in  a  preponderance  of  negative  values  for  mk.) 

Table  6.  Sample  Moment  Results  for  (75) 


Moment  k 

m^  for  N  *  100 

mk  for  N  *  1000 

mk  for  N  =  10000 

1 

1.07 

1.37 

0.52 

2 

1.64 

0.78 

0.65 

3 

1.08 

1.09 

0.42 

4 

1.49 

0.67 

0.89 

5 

0.95 

1.11 

0.30 

6 

1.21 

0.51 

0.94 

23 


TR  6843 

FIRST-ORDER  DISTRIBUTION  TEST 


Let  Ik  be  an  interval  of  the  line  segment  (0,1).  Then  let  cn  be  a 
counting  variable  given  by 


1  if  x  e  Ik 

c  =  S  >  f or  1  <  n  <  N 

'  0  otherwise 


(76) 


Then  an  estimate  of  the  probability  Pkt  that  xn  e  Ik,  is  given  by  the 
quantity 


i  4 

qk  =  N  nTlC" 


(77) 


This  random  variable  has  mean  and  variance 


E(qk)  =  Pk,  Var(qk)  =  Pk(l-Pk)/N 


(78) 


Thus  the  random  variable 


#k '  [yi-v'"]* 

I 


(79) 


should  lie  most  often  in  the  region  ±  2  if  xp  truly  does  have  probability 
Pk  of  lying  in  Ik.  In  table  7  is  displayed  the  results  of  sample  runs  for 
the  case  where  interval  (0,1)  is  broken  into  10  equal  parts;  that  is. 


(80) 


Once  again,  there  is  no  reason  to  reject  the  hypothesis  that  the  random  number 
generator  is  uniformly  distributed  over  (0,1). 
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Sample  Probability  Results  for  (79) 


Interval  k 


vk  for  N  =  1000 


v^  for  N  =  10000 


vk  for  N  =  100000 
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\  SUMMARY 

Several  methods  of  generating  random  variables  with  specified  cumulative 
distribution  functions  have  been  presented  and  evaluated  in  terms  of  their 
time  of  execution  and  efficiency  of  generation.  They  include  nonlinear 
distortion  of  a  single  (uniformly  distributed)  random  variable  or  through 
combinations  of  simply  generated  random  variables.  The  former  approach 
requires  the  ability  to  realize  the  inverse ^function  to  the  desired  cumulative 
distribution  function,  whereas  the  latter  approach  is  very  fruitful  if  the 
characteristic  function  corresponding  to  the  specified  cumulative  distriDution 
function  can  be  broken  down  into  a  product  of  simpler  components.  Which 
approach  to  adopt  depends  on  the  particular  example  and  the  exact  way  that  any 
parameters  enter  into  the  cumulative  distribution  function.  Of  course,  the 
inverse  to  a  cumulative  distribution  function  can  always  be  evaluated 
numerically  and  used  as  a  table  look-up,  in  order  to  generate  transformed 
random  variables;  however,  this  numerical  procedure  would  have  to  be  repeated 
if  any  parameter  of  the  cumulative  distribution  function  were  desired  changed. 

\ 

The  ability  to  generate  uniformly  distributed  random  variables  is  a  key 
component  of  this  procedure.  Several  statistical  tests  on  the  Hewlett-Packard 
9845  random  number  generator  have  confirmed  it  to  have  a  flat  probability 
density  function  and  independent  samples. 

These  techniques  are  useful  for  generating  sets  of  random  variables  of 
size  N  with  a  specified  cumulative  distribution  function,  and  then  plotting 
the  sample  cumulative  distribution  function* as  described  in  [6],  in  order  to 
ascertain  the  amount  of  fluctuation  which  is  typical  for  that  set  size  N  and 
in  different  regions  of  probability.  Then  by  comparing  the  amount  of 
fluctuation  of  a  measured  data  set  (of  unknown  statistics)  with  typical  sample 
cumulative  distribution  functions  of  the  same  size  and  with  a  known  specified 
cumulative  distribution  function,  decisions  on  acceptance  or  rejection  of  a 
hypothesized  cumulative  distribution  function  can  be  made  witn  confidence. 

See  [6,  figures  6  and  7]  for  illustrations  of  this  procedure. 
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APPENDIX  A.  SOME  INVERSE  FUNCTION  RELATIONS 
In  this  appendix,  x  and  y  are  real  (non  random)  variables.  Let  g(x)  be  a 


monotonic  function  for  x  in  a  given  range,  and  let 

y  -  9(x)  •  ( A— 1 ) 

The  inverse  relation  to  ( A— 1 )  is,  for  y  in  the  appropriate  range, 

x  =  9(y)  .  (A-2) 

Now  suppose  that  we  cascade  nonlinearities:  let 

z  =  h(y)  =  h(g(x) )  s  f(x)  ,  (A-3) 

to  yield  overall  nonlinearity  f.  Then  the  outer  equality  in  (A-3)  yields 

x  -  f(z)  ,  (A-4) 

whereas  the  first  and  third  terms  in  (A-3)  yield 

h(z)  «  g(x),  and  g(fi(z))=x  .  (A-5) 


Combining  (A-4)  and  (A-5),  there  follows  for  the  inverse  function  of  the 
cascade,  (A-3), 

f(z)  =  g(h(z) )  ,  (A-6) 

in  terms  of  the  inverses  of  the  individual  transformations. 

If  we  combine  (A-l)  and  (A-2),  we  get 

y  =  g(3(y))  •  (a-7) 

Taking  a  derivative  with  respect  to  y,  we  find 
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i  =  g'(g(y) )  g‘(y)  ; 

?'(y)  =  iW  •  1*'8) 

That  is,  the  derivative  of  the  inverse,  g,  can  be  found  in  terms  of  the 
inverse  function  and  the  derivative  of  the  original  function  g. 

Suppose  that  y  is  given  in  terms  of  x  via  transformation 

y  =  g(?(x) )  ,  (a-9) 

but  the  inverse  function  f  is  impossible  or  very  difficult  to  obtain  from 
given  function  f.  A  simple  way  of  evaluating  y  vs  x,  then,  is  parametrically 
by  letting 


t  *  f { x ) ,  to  get  x  =  f(t),  y  =  g(t) 


(A-10 ) 


Now,  as  t  is  varied,  f  and  g  can  be  evaluated  to  determine  y  vs  x.  This 
transformation  also  allows  evaluation  of  an  integral  involving  an  inverse 
function: 


b  tb 

£  dx  g(f(x) )  «  §  dt  f'(t)  g(t)  ,  (A-ll) 

a  t 


where 


*3  =?(*)»  tb  =  ?(b)  .  (A-12) 

Another  use  of  inverse  functions  in  integral  evaluation  is  afforded  by  the 
example 


I 


b 

$  dx  g(x)  w( x ) 
a 


(A— 13) 
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where  function  w(x)  need  not  be  monotonic  or  possess  an  inverse.  g(x)  is 
assumed  monotonic  in  (a,b).  Letting 


y  =  g(x),  X  =  g(y) ,  dx  =  dy  g'(y) 


there  follows,  for  the  integral  in  terms  of  inverse  function  <J, 


(A-14) 


g(b) 

i  =  5  dy  y  3'(y)  w(3(y)) 
g(a) 

Integrating  by  parts,  using  identification 


(A— 15) 


U  =  y  w(3(y)),  V  =  g(y) 


(A-16) 


we  get  the  alternative  form 


I  =  b  g(b)  w(b)  -  a  g(a)  w(a)  -  J  dy  3(y)  3-  {y  w(g(y))}  .  (A-17) 

g(a) 

The  special  case  of  w(x)  =  1  in  ( A— 1 3)  and  (A-17),  namely 


y  v u ) 

^  dx  g(x)  =  b  g(b)  -  a  g(a)  -  J  dy  g(y) 
a  g(a) 


(A— 18) 


has  the  geometrical  interpretation  in  figure  A-l.  Equation  (A— 18 )  is  the 
statement  that  Aj  +  ^  +  A3  =  total  area  b  g(b). 

As  an  example  of  the  application  of  (A-13)  and  (A— 15) ,  consider 


g(x)  =  arc  sin(x),  g(y)  =  sin  y 


(A-iy) 


There  follows 


b  arc  sin(b) 

^  dx  arc  sin(x)  w(x)  =  J  dy  y  cos(y)  w(sin  y) .  ( A— 20 ) 

a  arc  sin(a) 
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So  if  w  is  a  polynomial,  the  integral  on  the  right  side  of  ( A— 20 )  can  be 
evaluated  in  closed  form. 


Figure  A-l.  Geometrical  Interpretation  of  (A— 18) 
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APPENDI.'.  B.  A  MORE  GENERAL  DISTORTION 

In  (21),  a  random  variable  r  with  a  Rayleigh  probability  density  function 
and  a  random  variable  e  with  a  uniform  probability  density  function  were 
generated.  These  were  then  used  in  nonlinear  transformation  (23)  to  generate 
two  independent  Gaussian  random  variables  y^  and  y2.  Here  we  will 
generalize  the  probability  density  function  pr  0f  r  in  (22),  and  allow  pr 
to  be  arbitrary,  e  is  still  uniform  over  2 w . 

The  two  new  random  variables  generated  are  again  as  in  (23): 

yi  =  r  cos  e,  y2  =  r  sin  e  .  (B-I) 

Because  of  the  uniform  distribution  of  e  over  2n,  the  joint  probability 
density  function  of  y^  and  y2  is  of  the  symmetric  form 


y2)  =  +  y2  )  for  all  yv  y2 


( B— 2 ) 


To  determine  h,  observe  that,  for  t  >  0, 


P  *  Probj^  *  yf<  tj  .  j"jdyj  dy2  p(yj.  y2) 

-  dyl  h  (ft  +  yjQ  =  J  dp  p  d^  h(p)  =  2ir  ^ dp  p  h(p)  ,  { B— 3) 


where  Ct  iS  a  circle  of  radius  t  centered  at  the  origin.  But  also,  from 
(B-l) , 


P  =  Prob  fr  <  t]  =  J*du  p  (u)  .  (b-4) 

1  o 

Equating  (B-3)  and  (B-4)  and  taking  a  derivative  with  respect  to  t,  there 
follows 


Pr(t) 

h(t)  =  27E —  for  t  >  0 


(8-5) 
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Reference  to  (B-2)  then  yields  for  the  joint  probability  density  function  of 

yl.  y2> 


p(yx »  y2)  = 


for  all  yp  y2 


(B- 


yi  and  y2  are  statistically  dependent  in  general. 
EXAMPLE  1. 


Our  first  case  is  the  probability  density  function  in  (34),  for  random 
variable  r  as  generated  by  (33): 


Pr(u) 


u2N  1  exp(-u2/2) 
2N-1(N-1) J 


for 


u  >  0 


Substitution  in  (B-6)  yields  the  joint  probability  density  function 


(B- 


p(yp  y2)  -  [2ir(N-l).]_1 


for  all  yp  y2  . 


(B- 


The  special  case  of  N=1  reduces  to  the  pair  of  Gaussian  random  variables 
already  considered  in  (2 1 )— (24) .  For  N=2,  ( B— 7 )  and  ( B— 8)  yield 


pr(u)  -  \  u3  exp(-u2/2) 

2  2 

yf  +  y\ 

p(yi»  y2)  =  — ^ —  exp 

EXAMPLE  2. 


(B- 


Here  random  variable  r  is  generated  according  to  (32),  and  has  the 
probability  density  function  given  by  (29): 
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Substitution  of  (B-10)  in  { B— 6 )  yields 


N  _  x 

P(yr  y2)  =  [2*  (N-l).']-1  ^  exp(^fT^)  for  all  y^ 

Special  cases  for  N  equal  to  1  and  2  are  respectively 


and 


p(yr  y2)  =  (2ir)  1 


p(yi»  y2)  =  (2ir)_1 


y2.  ( B— l l ) 


(8-12) 


( B— 1 3 ) 
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APPENDIX  C.  A  TEST  FOR  WHITENESS  OF  A  SEQUENCE 


Suppose  data 


points  fxk) 


K-l 

0 


are  avai lable. 


estimate  at  delay  n  according  to* 


We  can  form  a  correlation 


5n'^xkxk-r>  fora""  •  (C-D 

where  summations  without  limits  are  over  the  range  of  nonzero  summands. 

Rn  is  nonzero  only  for  Jn|  <  K.  If  process  [xn]  were  white,  we  would 
expect  to  have 


| Rn |  <<  R0  for  all  n  ^  0  .  (C-2) 

Therefore  a  measure  of  nonwhiteness  is  afforded  by  the  ratio  E/ft*,  where  error 
measure  E  is  defined  as  the  sum  of  squares  for  n  ^  0, 


(C— 3 ) 


However,  (C-3)  is  very  time-consuming  to  calculate  via  ( C  — 1 ) ,  because  of  all 
the  multiplications  and  additions  required.  A  much  more  practical  evaluation 
of  (C-3)  is  afforded  by  the  following  procedure  (the  derivation  is  presented 
later  in  this  appendix). 


Define  an  M-point  OFT  of  the  K  data  points: 


K-l 

X  =  \  exp(-i2irkm/M)  for  0  <  m  <  M-l  .  (C-4) 

m  k=0 


Then  it  follows  that  if  M  >  2K-1, 


*  The  following  mathematical  development  is  very  similar  to  that  in  [7J. 
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and  our  whiteness  measure  can  be  expressed  as 


A  threshold  test  of  (C-6)  is  equivalent  to 


Threshold 


(< 


The  distribution  of  Q  could  be  computed  from  white-noise  simulations,  and  a 
threshold  value  selected  for  prescribed  error  probability.  By  Schwarz's 
inequality,  Q  >  1,  with  equality  realized  if  and  only  if  |XmJ^  is  constant 
for  all  m. 

Evaluation  of  Q  in  (C-7)  requires  one  M-point  DFT  of  the  data  {x^g-1, 
where  we  must  have 


M  >  2K  -  1 


0 


The  subsequent  calculations  in  (C-7)  are  quickly  conducted. 

An  alternative  interpretation  of  error  measure  E  in  (C— 3 )  is  very 
illuminating  and  lends  additional  credence  to  (C-7)  as  an  appropriate 
statistic.  If  we  define  spectral  estimate 

^  Rn  exp(-i2*mn/M)  for  all  m  ,  (( 

n 

then  we  find 

L  -  r  kJ2  for  0  <  m  <  M-l  ,  (C- 
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and  the  error  measure  E  in  (C— 3 )  becomes 


i  M-l  ,  A  M-l 


(C-ll) 


But  notice  that  if  we  define  the  sample  mean  of  the  set  of  spectral 
r*  1  M-l 


estimates 


M 


,  M-l 


(C— 12) 


then  the  sample  variance  becomes 

M-l 

2.  1  K  »2 

a  5  M  (Gm  “ 


ISS2  A 
M^)  ra  \M  ^0  " 


( C— 1 3) 


That  is,  error  measure  E  is  equal  to  the  sample  variance  of  the  set  of 
spectral  estimates.  This  latter  quantity  is  a  very  natural  measure  for 
deducing  whether  a  sequence  is  white,  since  a ^  would  be  expected  to  be 
smaller  for  a  truly  white  process. 

This  also  suggests  that  a  meaningful  measure,  for  determining  whether  a 
sequence  is  white  over  a  limited  band  of  the  zero-to-Nyquist  range,  is  the 
sample  variance  of  the  spectral  estimates  over  that  particular  band  in 
question. 

For  a  real  sequence  the  sums  in  (C-7)  need  only  be  conducted  over 

half  of  the  range,  by  virtue  of  the  conjugate  symmetry: 


XM-m  *  Xm  ^or  1  <  m  <.  M-l  if  (xjJ  real 


(C-14) 
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r»  !■  i- -r"  t  -a- — - — 

.  -  - 

r 

Thus  then  (for  M  even) 

t 

r 

M  -  1 

lxm|4  *  X0  *  4  *  2  Kl4  1f  l\]  real- 

(C— 15) 

r 

4 

* 

1 

9 

DERIVATIONS 

Define  spectral  estimate 

r* 

■  i 

exp(-i2irtim/M)  for  all  m 

( C  —1 6 ) 

..  J 

.] 

n 

J 

. 

The  superscript  specifically  indicates  the  period  in  tne  subscript  variable. 

Substituting  (C— 1 )  in  (C— 16) ,  there  follows 

=  2exP(~i2l,nm/M)  K  -eE.  \  xjC_n* 

i 

« 

=  ^  ^xk  exp(-i2itkm/M)  xk  n  exp(-i2ir(k-n)m/M)j 

* 

i 

r; 

=  I  |x(M)|2 

K  Tm  I 

(C— 17) 

5 

where 

r  < 

xiM^  *  IE.  K  exp(-i2irkm/rt)  for  all  m 
m  k  k 

(C— 18) 

;j 

-  • 

Now  the  inverse  DFT  of  is 

• 

*■  - 

^M)=  M-jg  §iM)  exp(i2wmn/M)  a  ...  *  Rn41  *  Rn  +  Rn+M  ♦  .. 

.  (C-19) 

j 

A 

■  i 
.  * 

i 

i-i 

•  v 

for  all  n, 

which  is  an  infinitely-aliased  version  of  [Rn] .  However,  if 

—a 

• 

i  : 

M  >  2K-1, 

38 

(C-20) 

— i 

* 

»■  ■ 

.1  , 

F 

• 

y-v  v  v—  ■■■  "jw  ■'— «  T*'  -T^T^  11 


kw- 

K.* 


& 


■ 
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there  is  no  overlap  of  terms  in  (C-19).  Henceforth  we  presume  (C— 20 )  to  be 
true;  that  is,  the  size  of  the  transform  (C-18)  must  be  at  least  twice  as 
large  as  the  number  of  data  points.  Then  there  follows  from  (C-19), 


r  _  I  xj  £(M) 

n  -  M  ^  fim 


exp  (i2*nm/M)  for  |n|  <  K 


Therefore 


5  k 


M-l  M-l  M-l 


n=0  ra=0  p=0 


K-l 

lfinl 

I2-  ^ 

1  “  n:one  period 

n=-K+l 

of  length  M 

M-l 

g(M)  g(M) 

m  p 

exp(i2n(m-p)n/M)  = 

*2 

m=0 

Then  from  (C— 3 )  and  (C-19), 

M-l 


M-l 


'-lie- »§«•< 


•)'  • 


( C— 2 1 ) 


( C  —2  2 ) 


(C— 23) 


Although  (C-18)  is  defined  for  all  relative  sizes  of  M  and  K,  the  relation 
(C— 23 )  holds  only  if  (C-20)  holds. 
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